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Disclaimer 


None of the authors, nor any contributor to the JETSPIN code or its derivatives guarantee that the 
software and associated documentation is free from error. Neither do they accept responsibility for any loss 
or damage that results from its use. The responsibility for ensuring that the software is fit for purpose lies 
entirely with the user. 
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About JETSPIN 


The code is licensed under Open Software License v. 3.0 (OSL-3.0). The full text of the licence can be 
found on the website: http://opensource.Org/licenses/OSL-3.0 

If results obtained with this code are published, an appropriate citation would be: 

Marco Lauricella, Giuseppe Pontrelli, Ivan Coluzza, Dario Pisignano, Sauro Succi, JETSPIN: A specific- 
purpose open-source software for electrospinning simulations of nanofibers, Submitted, to Computer Physics 
Communications , 2015. 
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Chapter 1 

Introduction 


Scope of Chapter 

This chapter describes the design and directory structure of JETSPIN. 

1.1 The JETSPIN Package 

In the recent years, electrospun nanofibers have gained a considerable industrial interest due to many possible 
applications, such as tissue engineering, air and water filtration, drug delivery and regenerative medicine. In 
particular, the high surface-area ratio of the fibers offers an intriguing prospect for technological applications. 
As consequence, several studies were focused on the characterization and production of uni-dimensionally 
elongated nanostructures. A number of reviews mmmmm and books U0! concerning electrospinning 
have been published in the last decade. 

Typically, electrospun nanofibers are produced at laboratory scale by the uniaxial stretching of a jet, 
which is ejected at the nozzle from a charged polymer solution. The initial elongation of a jet can be 
produced by applying an externally electrostatic field between the spinneret and a conductive collector. 
Electrospinning involves mainly two sequential stages in the uniaxial elongation of the extruded polymer jet: 
an initial quasi steady stage, in which the electric field stretches the jet in a straight path away from the 
nozzle of the ejecting apparatus, and a second stage characterized by a bending instability induced from 
small perturbations, which misalign the jet from its axis of elongation |9j. These small disturbances may 
originate from mechanical vibrations at the nozzle or from hydrodynamic-aerodynamic perturbations within 
the experimental apparatus. Such a misalignment provides an electrostatic-driven bending instability before 
the jet reaches the conductive collector, where the fibers are finally deposited. As a consequence, the jet 
path length between the nozzle and the collector increases, and the stream cross-section undergoes a further 
decrease. The ultimate goal of electrospinning process is to minimize the radius of the collected fibers. By 
a simple argument of mass conservation, this is tantamount to maximizing the jet length by the time it 
reaches the collecting plane. By the same argument, it is therefore of interest to minimize the length of the 
initial stable jet region. Consequently, the bending instability is a desirable effect, as it produces a higher 
surface-area-to-volume ratio of the jet, which is transferred to the resulting nanofibers US- 

Computational models provide a useful tool to elucidate the physical phenomenon and provide information 
which might be used for the design of electrospinning experiments. Numerical simulations can be used to 
improve the capability of predicting the key processing parameters and exert a better control on the resulting 
nanofiber structure. In recent years, with a renewed interest in nanotechnology, electrospinning studies 
attracted the attention of many resarchers both from modeling, computational and experimental point of 
view HU EH E3J M ESI M- Although some author uses dissipative particle dynamics (DPD) mesoscale 
simulation method into electrospinning study HSj, most models usually treat the jet filament as a series of 
discrete elements (beads) obeying the equations of continuum mechanics [111 112] . Each bead is subject to 
different types of interactions, such as long-range Coulomb repulsion, viscoelastic drag, the external electric 
field. The main aim of such models is to explain the complexity of the resulting dynamics and provide the 
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set of parameters driving the optimal process. The effect of fast-oscillating loads on the bending instability 
have been explored in an extensive modeling and computational study m- 

JETSPIN delivers a FORTRAN code especially designed to simulate the electrospinning process in a 
variety of different conditions and experimental settings. This comprehensive platform is devised such that 
different cases and input variables can be described and simulated. The framework is developed to exploit 
several computational architectures, both serial and parallel. 

JETSPIN, as open-source software, can be used to carry out a systematic sensitivity analysis over a broad 
range of parameter values. The results of simulations provide valuable insight on the physics of the process 
and can be used to assess experimental procedures for an optimal design of the equipment and to control 
processing strategies of technological advanced nanofibers. 


1.2 Structure 


JETSPIN is written in free format FORTRAN90, and it consists of approximately 160 subroutines. The 
source exploits the modular approach provided by the programming language. All the variables having in 
common description of certain features or method are grouped in modules. The convention of explicit type 
declaration is adopted, and all the arguments passed in calling sequences of functions or subroutines have 
defined intent. We use the PRIVATE and PUBLIC accessibility attributes in order to decrease error-proneness 
in programming. 

The main routines have been gathered in the main, f90 file, which drives all the CPU-intensive computa¬ 
tions needed for the capabilities mentioned below. The variables describing the main features of nanofibers 
(position, velocity, etc.) are declared in the nanojet_mod.f90 file, which also contains the main subroutines 
for the memory management of the fundamental data of the simulated system. Since the size system is 
strictly time-dependent, JETSPIN exploits the dynamic array allocation features of FORTRAN90 to assign 
the necessary array dimensions. In particular, the size system is modified by the routines add_jetbead and 
erase_jetbead, while the decision of the main array size, declared as mxnpjet, is handled by the routine 
reallocate_jet. The sizes of various service bookkeeping arrays are handled within a parallel implementa¬ 


tion strategy, witch exploits few dedicated subroutines (see Sec 1.31. All the implemented time integrators 


are written in the integrator_mod.f90 file, which contains the routine driver_integrator to select the 
proper integrator, as indicated in the input file. All the terms of equations of motion for the implemented 


model (see Sec 2.1 ) are computed by routines located in the eom_mod.f90 file, which call other subroutines 
in the files coulomb_force_mod.f90, viscoelastic_force_mod. f90 and support_functions_mod. f90. A 
summarizing scheme of the main JETSPIN program in the main, f90 file has been sketched in Fig o 

The user can carry out simulations of nanofibers without a detailed understanding of the structure of 
JETSPIN code. All the parameters governing the system can be defined in the input file (see Sec |3.1[ ), 
which is read by routines located in the io_mod.f90 and parse_mod.f90 files. Instead, the user should be 
acquainted with the model described in Cap [2j The content of the output file is completely customizable 
by the input file as described in Sec |3.1[ and it can report different time-averaged observables computed by 
routines of the module statistic_mod (see Sec |3.2| ). The routines in the error_mod. f 90 file can display 
various warning or error banners on computer terminal, so that the user can easily correct the most common 
mistakes in the input file. 

JETSPIN is supplied as a single UNIX compressed (tarred and gzipped) directory with four sub-directories. 
All the source code files are contained in the source sub-directory. The examples sub-directory contains dif¬ 
ferent test cases that can help the user to edit new input files. The build sub-directory stores a UNIX 
makefile that assembles the executable versions of the code both in serial and parallel version with different 
compilers. Note that JETSPIN may be compiled on any UNIX platform. The makefile should be copied 
(and eventually modified) into the source sub-directory, where the code is compiled and linked. A list of 
targets for several common workstations and parallel computers can be used by the command "make target 
”, where target is one of the options reported in Tab |l.l| On Windows system we advice the user to compile 
JETSPIN under the command-line interface Cygwin [16] . Finally, the binary executable file can be run in 
the execute sub-directory. 
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target: 

gfortran 

gfortran-mpi 

cygwin 

cygwin-mpi 


intel 

intel-mpi 

intel-openmpi 

help 


meaning: 

compile in serial mode using the GFortran compiler. 

compile in parallel mode using the GFortran compiler and the Open Mpi library, 
compile in serial mode using the GFortran compiler under the command-line 
interface Cygwin for Windows. 

compile in parallel mode using the GFortran compiler and the Open Mpi library 
under the command-line interface Cygwin for Windows (note a precompiled 
package of the Open Mpi library is already available on Cygwin). 
compile in serial mode using the Intel compiler. 

compile in parallel mode using the Intel compiler and the Intel Mpi library, 
compile in parallel mode using the Intel compiler and the Open Mpi library, 
return the list of possible target choices 


Table 1.1: List of targets for several common workstations and parallel computers, which can be used by the 
command "make target". 


1.3 Parallelization 

The parallel infrastructure of JETSPIN incorporates the necessary data distribution and communication 
structures. The parallel strategy underlying JETSPIN is the Replicated Data (RD) scheme [15], where funda¬ 
mental data of the simulated system are reproduced on all processing nodes. In simulations of nanofibers, the 
fundamental data consist of position, velocity, and viscoelastic stress arrays of each bead in which the nanojet 
is discretised (see below Cap[2|. Further data defining mass and charge of each bead are also replicated. 
However, all auxiliary data are distributed in equal portion of data (as much as possible) for each processor. 
Despite being available other parallel strategies such as the Domain Decomposition [SO], our experience has 
shown that such volume of data is by no means prohibitive on current parallel computers. 

By the RD scheme, we adopt in JETSPIN a simple communication strategy of the communication be¬ 
tween nodes, which is handled by global summation routines. The module version_mod located in the 
parallel_version_mod.f90 file contains all the global communication routines which exploit the MPI 
(Message Passing Interface) library. Note that a FORTRAN90 compiler and an MPI implementation for 
the specific machine architecture are required in order to compile JETSPIN in parallel mode. An alternative 
version of the module version_mod is located in the serial_version_mod.f90 file, and it can be easily 
selected by appropriate flags in the makefile at compile time. By selecting this version, JETSPIN can also 
be run on serial computers without modification, even though the code has been designed to run on parallel 
computers. 

The size system is strictly time-dependent as mentioned in Sec |1.2| and, therefore, the memory of vari¬ 
ous service bookkeeping arrays is dynamically distributed over all the processing nodes. In particular, the 
bookkeeping array size, declared as mxchunk, is managed by the routine set_mxchunk. All the beads of 
the discretised nanofiber are assigned at every time step to a specific node by the routine set_chunk, and 
their temporary data are stored in bookkeeping arrays belonging to the assigned node. It is worth stressing 
that the communication latency makes the parallelization efficiency strictly dependent on the system size. 
Therefore, we only advice the use of JETSPIN in parallel mode whenever the user expects a system size with 
at least 50 beads for each node (further details in Ref [21]). 
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Figure 1.1: Structure of the main JETSPIN program. 
































Chapter 2 


JETSPIN Model and Algorithms 


Scope of Chapter 

This chapter describes the model and simulation algorithms incorporated into JETSPIN. 

2.1 Equations of Motion 

The model implemented in JETSPIN is an extension of the Lagrangian discrete model introduced by Reneker 
et al. EJ (in the following text referred as Ren model). The model provides a compromise of efficiency and 
accuracy by representing the filament as a series of n beads (jet beads) at mutual distance l connected by 
viscoelastic elements. The length l is typically larger than the cross-sectional radius of the filament, but 
smaller than the characteristic lengths of other observables of interest (e.g. curvature radius). Each i — th 
bead has mass m, and charge qi (not necessarily equal for all the beads). The nanofiber is modelled as a 
viscoelastic Maxwell fluid, so that the stress eq for the element connecting bead i with bead i +1 is given by 
the viscoelastic constitutive equation: 


diji G dli G . 

—— = - - - -<7j, (2-1) 

dt li dt fj, 

where U is the length of the element connecting bead i with bead * + 1, G is the elastic modulus, /i the 
viscosity of the fluid jet, and t the time. Given cq the cross-sectional radius of the filament at the bead i, the 
viscoleastic force f ve pulling bead i back to * — 1 and towards i + 1 is 

f); e , i — 'Kd^Oi • t.j T (2-2) 

where t, is the unit vector pointing bead i from bead i — 1. The surface tension force f S ( for the i — th bead 
is given by 


f st,i — ki ’ 7T 


f Cli + di -1 

V 2 


2 

Oi ‘ C i , 


(2.3) 


where a is the surface tension coefficient, ki is the local curvature, and Cj is the unit vector pointing the 
center of the local curvature from bead i. Note the force f s * is acting to restore the rectilinear shape of the 
bending part of the jet. 

In the electrospinning experimental configuration an intense external electric potential J>o is applied be¬ 
tween the spinneret and a conducting collector located at distance h from the injection point. As consequence, 
each i — th bead endures the electric force 


h, = ■ x, (2.4) 

h 

where x is the unit vector pointing the collector from the spinneret and assumed along the x axis. Note that 
in Ren model the electric field $ is assumed to be static in order to avoid the computationally expensive 
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integration of Poisson equation, whereas in reality $ is depending on the net charge of the nanofiber so as to 
maintain constant the potential at the electrodes. The latter issue was elegantly addressed by Kowalewski et 
al. [22], and its implementation in JETSPIN will be planned. 

The net Coulomb force f c acting on the i — th bead from all the other beads is given by 


where Rij = 


(Xi - Xj ) 2 + (yi - Uj ) 2 + [Zi - ZjY 


j =i 13 

1 1/2 


(2.5) 


is the unit vector pointing the i — th bead 


and u 

from j — th bead.. Although the Ren model provides a reasonable description for the spiral motion of 
the jet, the last term f c introduces mathematical inconsistencies due to the discretization of the fiber into 
point-charges. Indeed, the charge induces a field on the outer shell of the fiber and not on the center line 
(as in the implemented model). Different approaches were developed to overcome this issue which usually 
imply strong approximations [251 1101 [2¥. . Other strategies use a less crude approximation by accounting 
for the actual electrostatic form factors between two interacting sections of a charged fiber [25] or involving 
more sophisticated methods which exploit the tree-code hierarchical force calculation algorithm [22] ■ The 
implementation in JETSPIN of methods based on tree-code hierarchical force calculation algorithm will be 
considered in future. 

Although usually much smaller than the other driving forces, the body force due to the gravity is computed 
in the model by the usual expression 


f g ,i = nmg • x, (2.6) 

where g is the gravitational acceleration. 

As shown in experimental evidences m, the air drag force affects the dynamics of the nanofiber. The 
extended stochastic model developed in Refs mm was included in JETSPIN to reproduce the aerodynamic 
effects. In particular, the air drag is modelled by adding two force terms: a random term and a dissipative 
term. The dissipative air drag term is usually dependent on the geometry of the jet, which changes in time, 
and it combines longitudinal and lateral components. Based on experimental results[261. the longitudinal 
component of the air drag dissipative force term acting on bead i is 

f air,i = -0.657 TUiliPa (2a l /^ a ) _0 ' 81 • tj_i (2.7) 

where p a and v a are the air density and kinematic viscosity, respectively, tj_i is the unit vector pointing 
bead i from bead i — 1, and Vtan,i — 33 i ■ tj_i is the projection of the vector velocity of bead i along the unit 
vector tj_i. Further, it is possible to write f a ir,i as 


f air,i — Z 


0.905 1+0.19 I* 
V t.i ' D-1 


if we introduced the dissipative coefficient y, equal to 


( 2 . 8 ) 


li 


0 65^0.905 Ay 

rrii 



(2.9) 


where V); = 7raf k is the jet volume represented by the i — th bead. 

Following the expression introduced by A. Yarin[2!jj||)|, under a high-speed air drag the lateral component 
of the aerodynamic dissipative force related to the flow speed is given in the linear approximation (for 
small bending perturbations) by 


f— ki ' Pa^i 33 


di + di -1 


( 2 . 10 ) 


where hi is still the local curvature, and c, is the unit vector pointing the center of the local curvature from 
bead i. The combined action of such longitudinal and lateral components provide the dissipative force term 
acting on the i — th bead 
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f diss,i — f air,i T f■ 

Instead, the random force term for the i — th bead has the form 


( 2 . 11 ) 


fraud,i = \j2m 2 i D v ■ ( 2 . 12 ) 

where D v denotes a generic diffusion coefficient in velocity space (which is assumed constant and equal for 
all the beads), and r) i is a three-dimensional vector, whereof each component ?y is an independent stochastic 
process that is nowhere differentiable with < rj (H) rj (t 2 ) >= S (\t 2 — ti|), and < rj(t) >= 0. Note that for 
the sake of simplicity we assume rj = dui(t)/dt, where uj(t) is a Wiener process. 

The combined action of these forces governs the elongation of the jet according to the Newton’s equation 
providing a non-linear Langevin-like stochastic differential equation: 


dlJi —> 

TYli ^ f el,i T “1“ f ve,i T f st,i T f g,i T f diss,i f rand t 


where is the velocity of the i — th bead. The velocity Vi satisfies the kinematic relation: 


(2.13) 


drj 

dt 


= Vj 


(2.14) 


where is the position vector of the i — th bead, ?* = a^x + yiy + ZiZ. The three Eqs 2. 31 |2.13| and |2.14| form 
the set of equations of motion (EOM) governing the time evolution of system. It is worth stressing that Eq 
|2.13| recovers a deterministic EOM by imposing p a = 0 and D v = 0. 


2.2 Perturbations at the nozzle 


The spinneret nozzle is represented by a single mass-less point of charge q fixed at x = 0 (nozzle bead). 
Its charge q is assumed equal to the mean charge value of the jet beads. Such charged point can be also 
interpreted as a small portion of jet which is fixed to the nozzle. In JETSPIN it is possible to add small 
perturbations to the y n and z n coordinates of the nozzle bead in order to model fast mechanical oscillations 
of the spinneret }17). Given the initial position of the nozzle 

y n = A- cos (<p) 

(2.15a) 

Z n = A- sin (ip), 

(2.15b) 

the equations of motion for the nozzle bead are 


dy n 

- 7 ~ = -w-z n 
dt 

(2.16a) 

dz n 

ld =UJ - y - 

(2.16b) 


where A denotes the amplitude of the perturbation, while w and p are its frequency and initial phase, 
respectively. 


2.3 Jet insertion 

The jet insertion at the nozzle is modelled as follows. For the sake of simplicity, let us consider a simulation 
which starts with only two bodies: a single mass-less point fixed at x = 0 representing the spinneret nozzle, 
and a bead modelling an initial jet segment of mass nii and charge located at distance lstep from the 
nozzle along the x axis. Here, l ste p denotes the length step used to discretise the jet in a sequence of beads. 
The starting jet bead is assumed to have a cross-sectional radius do, defined as the radius of the filament 
at the nozzle before the stretching process. Applying the condition that the volume of the jet is conserved, 
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the relation ir af Z* = 7ragZ step is valid for any i — th bead. Furthermore, the starting jet bead has an initial 
velocity v s along the x axis equal to the bulk fluid velocity in the syringe needle. Once this traveling jet bead 
is a distance of 2 • l s t ep away from the nozzle, a new jet bead (third body) is placed at distance Z step from the 
nozzle along the straight line joining the two previous bodies. Let us now label i — 1 the farthest bead from 
the nozzle, and i the last inserted bead. The i — th bead is inserted with the initial velocity Vi = v s + Vd, 
where Vd denotes the dragging velocity computed as 


Vd = 


Vi-i - v s 


(2.17) 


Here, the dragging velocity should be interpreted as an extra term which accounts for the drag effect of 
the electrospun jet on the last inserted segment. Note that the actual dragging velocity definition was chosen 
in order to not alter the strain velocity term (1 /k-i) • (dli_i/dt) of Eq 2.1 before and after the bead insertion. 


2.4 Dynamic refinement 


In the original Ren model, each bead represent a constant quantity V c of jet volume, so that the relationship 
7 raj h = 7 Ta,Ql s te P = V c is valid for any i — th bead. Despite such relationship can be sometimes a useful 
advantage (as example, see Sec 2.51, the radius reduction eq of the electrospun nanofiber provides a significant 
increase of the discretization length Z,; value along the jet path (also greater than two order of magnitude). 
As consequence, it was observed in literature m that the jet discretization close the collector becomes 
rather coarse to model efficiently the nanofiber. To get around this problem a dynamic refinement can be 
used in JETSPIN in order to maintain the the length of the elements Z, below a prescribed characteristic 
length max Z max . Whenever a bead length Zj is larger than Z max , the nanofiber is refined by discretizing 
uniformly the jet at the length step l s t ep value given in input file. The discretization is performed by Akima 
spline interpolations nasi] of the main quantities describing the jet beads (positions, fiber radius, stress, 
velocities). In order to perform the interpolation we proceed as follows: First of all, we introduce the variable 
A € [0,1] to parametrize the jet path, where A = 0 identifies the nozzle, and A = 1 the jet at the collector. 
For each bead we compute its respective \ value by using the formula 


1 

A i = - - J2 lk ( 2 - 18 ) 

hath k=1 

where l pa th = h is the jet path length from the nozzle to the collector, and i — 1 and i = Nf, ea d s 

denote the jet bead at the nozzle and collector, respectively (i = 0 denotes the nozzle). The set of A,; values 
represent the mesh used to build the cubic spline. Given a generic quantity y, the data Vi = y (A,) are 
tabulated and used in order to compute the coefficients of the Akima spline , following the algorithm reported 
in Ref [31] . Then, we enforce a uniform parametrization of the nanofiber by imposing all the lengths of the 
elements equal to Z step . The new mesh A* is defined as 

AS = 0 A* = —^—, t = 1,2,..., N£ eads . (2.19) 

v beads 

where N£ eads = l pa th/lstep represents the number of jet beads in the new representation. Finally, the new 
values y{ A*) are computed for any i — th bead by the spline interpolation. The procedure is repeated for 
positions, fiber radius, stress and velocities of the jet beads in order to provide the respective quantities in 
the new mesh A*. It is worth to sressing that the volume V r = 7raf U represented by each bead is not anymore 
constant, and, consequently, Traf h ^ n agZ step . In the following text we will refer to this emended version of 
the Ren model as DyRen model. 


2.5 Dimensionless quantities 

In JETSPIN all the variables are automatically rescaled and stored in dimensionless units. In order to adopt 
a dimensionless form of the equations of motion, we use the dimensionless scaling procdure proposed by 
Reneker et al. m- We define a characteristic length 
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Ln — 


I q 2 

7 TCLqG 


= l 


step 


nalPy 


G 


( 2 . 20 ) 


where we write the charge q as Tra^lgtepPy, denoting py the electric volume charge density of the filament. 
Further, we divide the time t and the stress a by their respective characteristic scales reported in Tat j2.1| 
Applying the condition that the volume of the jet is conserved 7raf li = V) for any i — th bead, we rewrite the 
EOM as: 


do-j 

dt 


dr-i -> 


-£=~ n 

(2.21a) 

1 dli 


k dt a1 

(2.21b) 


—= =4> ■ x 
dt 


■ u ij — • tj + F vet i+iVi+i - ' r • 

j —i dtij H H+l 

iAi 


&i +1 


■*T 


y/Vi 

Vk 


Vi- 


i—1 


F • x — r /D-aos-r+o.rg z* 

r Q x - 1 * l i u t.an.i L *-l 


(2.21c) 


- A. 


hvl 


y/Vi 

Vk 


Vi-i 


y/2&v-rfi(t) 


where we used the dimensionless derived variables and groups defined in Tal |2.1| It is worth stressing that 
the viscoelastic and surface tension force terms are slightly different from the dimensionless form provided 
by Reneker et al. d. since we are considering the DyRen model presented in Sec |2.4[ In particular, we do 
not use the relation tux\ li = na^lstep since it is valid only for the Ren model (see discussion in Sec 2.41. 
Finally, the dimensionless EOM of the nozzle are given as 


dy A - - K . 

iT - IY S 

at 

d A _ ry - 

di - Ks ' y " 


(2.22a) 


with the dimensionless parameter K s defined in Tat|2.1| 


(2.22b) 
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Characteristic Scales 


Ln — l 


™ a Z Py 


step y q 

a 0 = G 


Dimensionless Quantities 


i = Jl 

1 L 0 

Oi 

<Ji = — 

o-o 

k^L o 
Vi 


Rjj 

- ’ll 

Vi = Vi — 

L 0 

_ _i2_ 

Vtan,i Vtan,i T 

-^0 


V = 


Ll 


Dimensionless Groups 


i>, = 

F — 

x ve.i — 


Qi^oF 2 
rriih Lq G 2 
LqF 2 


miG 

9F 2 


Qij 


Ai = 


F = 

9 LnG 2 

H = r 


K. = uj 


rrHLlG 2 

afj? 
rriiG 2 
F 


hz r "{i ' toL, 

A = 


0.905 

0 


Pa-Lp 

m, 


0.19 


Instep 


0„ = D, 


G 

step 


An 


4 

Ll 


Table 2.1: Definitions of the characteristic scales, dimensionless derived variables, and groups employed in 
the text. 


2.6 Integration schemes 

In order to integrate the homogeneous differential equations of motion we discretise time as L = tp + * At with 
i = 1,, nsteps, where nsteps denotes the number of sub-intervals. Note that in JETSPIN the time step At is 
automatically rescaled by the quantity r, in accordance with the mentioned dimensionless scaling convention. 
Three different integration schemes are available for deterministic simulations: the first-order accurate Euler 
scheme, the second-order accurate Heun scheme (sometimes called second-order accurate Runge-Kutta), and 
the fourth-order accurate Runge-Kutta scheme [35 ■ The user can select a specific scheme by using appropriate 
keys in the input file, as described in Sec |3.1| For the integration of stochastic simulations the explicit strong 
order scheme by Platen [MIESES] was implemented, whereof the order of strong convergence was evaluated 
equal to 1.5. This scheme avoids the use of derivatives by corresponding finite differences in the same way 
as Runge-Kutta schemes do for ODEs in a deterministic setting, and it is briefly summarized as follows. 

Let us consider a Brownian motion vector process X = {X t , t} of d-dimensional satisfying the stochastic 
differential equation 


— =a(t,X 1 ,...,X d ) +bdO (2.23) 

dt 

where a and b are vectors of d-dimensional usually called drift and diffusion vector coefficients, and f 1 (t) 
denotes a Wiener process. Denoted Y t k the approximation for the fc-th component of the vector X at time t, 
the integrator has the following form: 
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with the vector supporting values 
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Y± = Y t + a At ± b '/At 

\ (2.25) 

3>± = Y + ±b(Y + J '/At.. 

Here, Afi and AT indicate normally distributed random variables constructed from two independent 
N (0,1) standard Gaussian distributed random variables (Ui, U 2 ) by means of the following linear transfor¬ 
mation: 


AH = UW At 

AT = ^A t 3/2 (Ux + ^=U^j. 


(2.26) 



Chapter 3 

JETSPIN Data Files 


Scope of Chapter 

This chapter describes all the input and output files of JETSPIN. 


3.1 Input file 

In order to run JETSPIN simulations an input file has to be prepared, which is a free-format and case- 
insensitive. The input file has to be named input.dat, and it contains the selection of the model system, 
integration scheme directives, specification of various parameters for the model, and output directives. The 
input file does not requires a specific order of key directives, and it is read by the input parsing module. Every 
line is treated as a command sentence (record). Records beginning with the symbol # (commented) and 
blank lines are not processed, and may be added to aid readability. Each record is read in words (directives 
and additional keywords and numbers), which are recognized as such by separation by one or more space 
characters. 

As in the example given in input test, the last record is a finish directive, which marks the end of 


the input data. Before the finish directive, a wide list of directives may be inserted (see Tab 3.1). The 
key systype should be used to set the nanofiber model. In JETSPIN two models are available: 1) the 
one dimensional model similar to the model of Cap [2] but assuming the nanofiber to be straight along the 
x axis, and, therefore, neglecting the surface tension force f s *; 2) the three dimensional model described 
in Cap [2] Internally these options are handled by the integer variable systype, which assumes the values 


explained in Tab 3.1 Further details of the 1-D model can be found in Refs 133 El ■ A series of variables 
are mandatory and have to be defined. For example, timestep, final time, initial length, etc. (see 


underlined directives in Tab 3.11. A missed definition of any mandatory variable will call an error banner 


on the terminal. The user can use two possible ways to define the initial geometry of a nanofiber both given 
the mandatory directive initial length: 1) the discretization step length of nanofiber is specified by the 
directive resolution, which causes automatically the setting of the jet segments number (the number of 
segment in which the jet is discretised); 2) the jet segments number is declared by the directive points, while 
the value of the discretization step length is automatically set by the program (as in Example Test Case 3 
in Sec 4.3). The directive cutoff indicates the length of the upper and lower proximal jet sections, which 


interact via Coulomb force on any bead. It is worth stressing that the length value is set at the nozzle, so 
that the effective cutoff increases along the simulation as the nanofiber is stretched. 

The user should pay special attention in choosing the time integration step given by the directive 
timestep, whose detailed considerations are provided in Ref |2Tj. The reader is referred to Tab 3.1 for 
a complete listing of all directives defining the electrospinning setup parameters. Not all these quantities 
are mandatory, but the user is informed that whenever a quantity is missed, it is usually assumed equal 
to zero by default (exceptions are stressed in Tab 3.11. Note all the quantities have to be expressed in 


centimeter-gram-second unit system (e.g. charge in statcoulomb, electric potential in statvolt, etc.). 
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directive: 

meaning: 


airdrag yes 

active the inclusion of the airdrag force in the model 
(note that if yes the directive “integrator 4” is mandatory) 

airdrag airdensity / 
airdrag airviscosity / 
airdrag airvelocity / 
airdrag amplitude / 

mass density of the air surrounding the nanofiber (default: 0 ) 
kinematic viscosity of the air surrounding the nanofiber (default: 1 ) 
velocity of the air with respect to the z axis (default: 0 ) 
amplitude of the randomic force rescaled by 7 (default: 1 ) 

(note that the dissipation coefficient 7 is automatically 
computed by JETSPIN using the input data) 

collector distance / 

distance of collector along the x axis from the nozzle 
(note the nozzle is assumed in the origin) 

cutoff / 

length of the proximal jet sections interacting by Coulomb 
force (default: equal to the collector distance) 

density mass / 
density charge / 
dynam refin yes 
dynam refin threshold / 
dynam refin every / 
dynam refin start / 

mass density of the nanofiber 

electric volume charge density of the nanofiber 

active the dynamic refinement during the simultion 

set the threshold beyond which the refinement is applied 

the refinement is applied every / seconds 

the refinement is applied only if the simulation time is greater 
than or equal to / seconds 

elastic modulus / 
external potential / 
final time / 
finish 

gravity yes 
initial length / 
integrator i 

elastic modulus of the nanofiber 

electric potential between the nozzle and collector 

set the end time of the simulation 

close the input file (last data record) 

active the inclusion of the gravity force in the model 

length of the nanofiber at the initial time 

set the integration scheme. The integer can be T’ for 

Euler, ’2’ for Heun, ’3’ for Runge-Kutta, and ’4’ 

for the Platen scheme (only for airdrag force activated). 

inserting yes 
nozzle cross / 
nozzle stress / 
nozzle velocity / 
perturb yes 
perturb freq / 
perturb ampl / 
points i 
print list Si . .. 

inject new beads as explained in Subsec 
cross section radius of the nanofiber at the nozzle 
viscoelastic stress of the nanofiber at the nozzle 
bulk fluid velocity of the nanofiber at the nozzle 
active the periodic perturbation at the nozzle 
frequency of the perturbation at the nozzle 
amplitude of the perturbation at the nozzle 
number of segment in which the jet is discretised 
print on terminal statistical data as indicated by symbolic 
strings (see Tab |3.2| for detailed informations) 

print time / 
print xyz / 

print data on terminal and output file every / seconds 
print the trajectory in XYZ file format in centimeters 
every / seconds 

print xyz frame / 

print the jet geometry in a single XYZ file format in centimeters 
every / seconds (default name of files: frame%06d.xyz) 

print xyz maxnum i 

print the XYZ file format with i beads starting from the 
nearest bead to the collector (default: 100 ) 

print xyz rescale / 
printstat list si ... 

print the XYZ file format data rescaled by / 
print on output file statistical data as indicated by symbolic 
strings (see Tab 3.2 for detailed informations) 

removing yes 
resolution / 
seed i 

surface tension / 
system i 

remove beads at collector 
discretization step length of nanofiber 
seed control to the internal random number generator 
surface tension coefficient of the nanofiber 
set the nanofiber model. The integer can be equal to 
’1’ for select the ld-model or ’3’ for the 3d-model 

timestep / 
viscosity / 

set the time step for the integration scheme 
viscosity of the nanofiber 


Table 3.1: Here, we report the list of directives available in JETSPIN. Note i, /, and s denote an integer 
number, a floating point number, and a string, respectively. The underlined directives are mandatory. The 
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3.2 Output files 


A series of specific directives causes the writing of output files (see Tab |3.1| ). In JETSPIN three output files 
can be written: 1) the file statdat.dat containing time-dependent statistical data of simulated process; 2) 
the file traj .xyz reporting the nanofiber trajectory in XYZ format file; 3) the files frame"/ 0 06d. xyz reporting 
the nanofiber geometries in splitted XYZ format files (a single file for each frame). 

Various statistical data can be written on the file statdat. dat , and the user can select them in input using 
the directive printstat list followed by appropriate symbolic strings, whose the list with corresponding 
meanings is reported in Tab |3.2| The selected observables will be printed at time interval on the same line 
following the order specified in input file. In the same way, a list of statistical data can be printed on computer 


terminal using the directive print list followed by the symbolic strings of Tab 3.2 


The file traj .xyz is written as a continuous series of XYZ format frames taken at time interval, so that 
can be read by suitable visualization programs (e.g. VMD-Visual Molecular Dynamics [33], UCSF Chimera 
|35J, etc.) to generate animations. The number of elements contained in the file is kept constant equal to 
the value specified by the directive print xyz maxnum, since few programs (e.g. VMD) do not manage a 
variable number of elements along the simulations. If the actual number of beads is lower than the given 
constant maxnum, the extra elements are printed in the origin point. 
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keys: 

meaning: 

t 

unsealed time 

ts 

scaled time 

x, y, z 

unsealed coordinates of the farest bead 

xs, ys, zs 

scaled coordinates of the farest bead 

st 

unsealed stress of the farest bead 

sts 

scaled stress of the farest bead 

vx, vy, vz 

unsealed velocities of the farest bead 

vxs, vys, vzs 

scaled velocities of the farest bead 

yz 

unsealed normal distance from the x axis of the farest bead 

yzs 

scaled normal distance from the x axis of the farest bead 

mass 

last inserted mass at the nozzle 

q 

last inserted charge at the nozzle 

cpu 

time for every print interval 

cpur 

remaining time to the end 

cpue 

elapsed time 

n 

number of beads used to discretise the jet 

f 

index of the first bead 

i 

index of the last bead 

curn 

current at the nozzle 

cure 

current at the collector 

vn 

velocity modulus of jet at the nozzle 

VC 

velocity modulus of jet at the collector 

SVC 

strain velocity at the collector 

mfn 

mass flux at the nozzle 

mfc 

mass flux at the collector 

rc 

radius of jet at the collector 

rrr 

radius reduction ratio of jet 

ip 

length path of jet 

rip 

length path of jet divided by the collector distance 


Table 3.2: In JETSPIN a series of instantaneous and statistical data are available to be printed by selecting 
the appropriate key. Here, the list of symbolic string keys is reported with their corresponding meanings. 
Note by scaled we mean that the observable was rescaled by the characteristic values provided in Tal |2.1| 
The underlined keys correspond to data which are averaged on the time interval given by the directive print 
time in input file. Note by farest bead we mean the farest bead from the nozzle (with greatest x value). Note 
all the quantities are expressed in centimeter-gram-second unit system, excepted the dimensionless scaled 
observables. 








Chapter 4 

Example Simulations 


Scope of Chapter 

This chapter describes the standard test cases for JETSPIN, the input and output files for which are in the 
examples sub-directory. 


4.1 Test Case 1: 1-D case 

This input file (located in the input-1 sub-directory) provides the dimensionless parameter values Q = 12, 
V = 2 and F ve = 12, which have been already used as reference case in Refs HUES]. 


4.2 Test Case 2: 1-D case with bead insertion at the nozzle acti¬ 
vated 

This input file (located in the input-2 sub-directory) provides the dimensionless parameter values Q = 12, 
V = 2 and F ve = 12, as in the previous test case, but here we have activated the injection of new beads by 
the directive inserting yes in the input file. 


4.3 Test Case 3: Three-D simulations of a process leading to poly¬ 
mer nanofibers 

The electrospinning of polyvinylpyrrolidone (PVP) nanofibers is a prototypical process, which has been 
largely investigated in literature. ED ED Hi] By the input located in the input-3 sub-directory, we simulate the 
electrospinning process of PVP solutions. Then, the theoretical results predicted by the models are compared 
with the aforementioned experimental data. In particular, we reproduce an experiment in which a solution 
of PVP (molecular weight = 1300 kDa) is prepared by a mixture of ethanol and water (17:3 v:v), at a 
concentration of about 2.5 wt%. The applied voltage is in a range around 10 kV, and the collector is placed 
at distance 16 cm from the nozzle, which has radius 250 micron (further details are provided in Ref. [ 101 ST] 
). As rheological properties of such system we consider the zero-shear viscosity po = 0.2g/(cm • s) [HI. 142]. 
the elastic modulus G = 5 • 10 4 g/(cm • s 2 ) [13], and the surface tension a = 21.1 g/s 2 [H]. We use for 
the simulation a viscosity value p which is two order of magnitude larger than the zero-shear viscosity po 
reported, since the strong longitudinal flows we are dealing with can lead to an increase of the extensional 
viscosity from po, as already observed in literature. Emn] The mass density is equal to 0.84 g/cm 3 , while 
the charge density was estimated by experimental observations of the current measured at the nozzle. For 


convenience, all the simulation parameters are summarized in Tat 4.1 
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p 

(kg/m 3 ) 

Pq 

(C/L) 

ao 

(cm) 

v s 

(cm/s) 

a 

(mN/m) 

P 

(Pa-s) 

G 

(Pa) 

Vo 

(kV) 

UJ 

(s- 1 ) 

A 

(cm) 

840 

2.8 • 10~ 7 

5•10 -3 

0.28 

21.1 

2.0 

50000 

9.0 

10 4 

10 -3 


Table 4.1: Simulation parameters for the simulation of PVP nanofibers. The headings used are as follows: 
p\ density, p q \ charge density, ao : fiber radius at the nozzle, v s : bulk fluid velocity in the syringe needle, a: 
surface tension, p: viscosity, G : elastic modulus, Vo : applied voltage bias, u>: frequency of perturbation, 
A : amplitude of perturbation. The bulk fluid velocity v s was estimated considering that the solution was 
pumped at constant flow rate of 2 mL/h in a needle of radius 250 micron. 
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